A Drosophila systems model of withdrawal from chronic pentylenetetrazole relevant in post - epileptogenesis

Background Rodent kindling induced by pentylenetetrazole (PTZ) is an established model of epileptogenesis and antiepileptic drug (AED) testing. Recently, a Drosophila systems model has been described in which chronic PTZ causes a decreased climbing speed in adult males on 7 day. Some AEDs ameliorate development of this locomotor deficit. Time-series of microarray expression profiles of heads of flies treated with PTZ has been found to resemble transcriptomic alterations associated with epileptogenesis. In the fly model, withdrawal from seven day long PTZ treatment causes an increased climbing speed on 7 consequent day. Here, we present a systems model of the post-PTZ withdrawal regime. Results Unlike AED-untreated individuals, flies treated with any of the five AEDs after PTZ discontinuation exhibited normal climbing speed on 7 day, i.e., 14 day from the beginning of PTZ treatment. Time-series of microarray expression profiles of fly heads comparing control PTZand AED-untreated, and AED-untreated post PTZ withdrawal groups showed differentially expressed genes throughout. These genes enriched gene ontology (GO) molecular functions including transcription regulator and GTPase regulator activities. Interestingly, expression profiles of fly heads comparing control PTZand AED-untreated, and AED-treated post PTZ withdrawal groups showed neutralization of transcription regulator and GTPase regulator activities by the AEDs. Further transcriptomic analysis based on overinteraction in protein interactome and enrichment of miRNA targets implicated axon guidance and neuronal remodeling related perturbations in the fly model. Conclusions Differential expression of genes belonging to transcription regulator and GTPase regulator activities have previously been reported in post-epileptogenesis in established rodent models. Also, axon guidance and neuronal remodeling related alterations have been implicated in epilepsy. The Drosophila model thus provides a unique opportunity to dissect long-term plasticity relevant in epileptogenesis at cellular and molecular levels. Besides, the model also offers an excellent system to efficiently screen agents with potential therapeutic activity. Background Epileptogenesis, the process of epilepsy development, is a network problem that may involve molecular, structural, and functional alterations in the brain. It is poorly understood in cellular and molecular terms [1, 2]. Overlapping pathophysiological mechanisms are considered to underlie epilepsy and various other neurological and psychiatric conditions [3, 4]. Besides epilepsy, AEDs are also used in treating these other nonepileptic conditions [3, 5, 6]. Rodent kindling induced by PTZ is one of the established models of epileptogenesis and AED testing [7]. Repeated injection of a subconvulsant dose of the GABA antagonist over several weeks results in partially and fully kindled animals with the latter group showing clonic-tonic seizures; once kindled, the state of behavioral hyperexcitability persists for up to several weeks after discontinuing PTZ [8, 9]. Understanding epileptogenesis at systems level may facilitate development of novel antiepileptogenic, disease-modifying, and neuroprotective agents [10]. The inherent complexity of mammalian brain however precludes systems modeling of epileptogenesis in rodents [10]. Recently, a Drosophila systems model of PTZ induced locomotor plasticity responsive to AEDs has been empirically developed [11]. Altered locomotor activity is one of the kindling and post-kindling behavior in rats and mice [12, 13]. Also, locomotor behavior of Drosophila adults has been used previously to model neuropsychiatric conditions [14]. Given this, the fly PTZ model is considered useful in understanding processes relevant in epileptogenesis and in testing AEDs [11]. In one aspect of this model, seven days of PTZ treatment causes a decrease in climbing speed of adult male flies. Concomitant treatment with the AEDs sodium valproate (NaVP) or levetiracetam (LEV), not ethosuximide (ETH), gabapentin (GBP) or vigabatrin (VGB), suppresses appearance of decreased climbing speed on 7 day [11]. Time-series of microarray expression profiles has revealed wide-spread transcriptomic alterations in the head of PTZ treated flies [11]. The convulsant drug mainly downregulated gene expression, with affected genes overrepresenting various GO biological processes like transcription, neuron morphogenesis during differentiation, synaptic transmission, regulation of neurotransmitter levels, neurogenesis, axonogenesis, protein modification, axon guidance, actin filament organization etc. [11]. Proteomic interactome based analysis has provided directionality to these events [11]. Further, pathway overrepresentation analysis has revealed enrichment of Wnt signaling and other associated pathways in genes downregulated by PTZ. Most importantly, mining of available genetic, transcriptomic or proteomic data pertaining to Drosophila seizure mutants, established rodent epiletogenesis models or human epileptic patients has shown overrepresentation of epilepsy associated genes in the PTZ regulated set in the fly model [11]. The post-PTZ withdrawal regime of the fly model has not been characterized. Here, we describe the behavioral and transcriptomal pharmacology of this regime. Our results suggest that postPTZ withdrawal regime models post-epileptogenesis. Methods Fly behavior and drug treatment Previously described [11] methods were used. In brief, 3-4 days old unmated males were first treated with 8 mg/ml of PTZ for 7 days and then shifted to drug free media for 7 days for climbing speed measurements, and for 1, 3 or 7 days for time-series microarray gene expression analysis. For the last time-point, flies were shifted to drug free media for three days and then again to vials containing fresh batch of drug free media. Flies treated with normal food (NF) were used as controls. The above three time-points correspond to 8, 10 and 14 days from the beginning of PTZ treatment, in that order. In AED treatment after withdrawing PTZ, flies were treated for first 3 days with an AED, i.e., up to 10 day, and the next four days with NF. Climbing speed was measured at the end, i.e., on 14 day of the beginning of PTZ treatment. Climbing speed was measured at the end, i.e., on 14 day of the beginning of PTZ treatment. Control flies (NF) were never exposed to any drug. Whereas climbing speed was measured only on 14 day, the expression profiles were generated for 10 and 14 day flies. Final concentration of PTZ, ETH, NaVP, GBP, VGB (all from Sigma-Aldrich) and LEV (Levesam 500, Nicholas Piramal) in the fly medium was 8, 3.48, 0.33, 16, 24 and 5 mg/ml, in that order. Climbing speed was measured using an indigenously developed semi-manual method that was validated by Dynamic Image Analysis System (DIAS v. 3.4.2, Soll Technologies) [11]. Student’s t-test, with Bonferroni corrected p-values, was used in the statistical analysis of climbing speed. Microarray profiling and data analysis Previously described [11] methods were performed. Total cellular RNA was isolated from fly heads belonging to four biological replicates. Microarray -cDNA Synthesis Kit, Target Purification Kit, and -RNA Target Synthesis Kit (Roche) were used to generate labeled antisense RNA. Starting with 10 μg of total cellular RNA, Eberwine method (kits from Roche) was used to generate cDNA and thereafter Cy and Cy (Amersham) labeled antisense RNA. The Cy and Cy labeled aRNAs (control and treated) were pooled together and precipitated, washed, air-dried, and dissolved in 18MΩ RNAase free water. A total of 52 microarrays (12Kv1, 44 nos.; 14Kv1, 8 nos.; CDMC) were hybridized, four each for 8 day and 7 day of PTZ treatment. Out of four, two slides were dye-swaps. Slides were scanned at 10 μm resolution using GenePix 4000A Microarray Scanner (Molecular Devices) and the images preprocessed and quantified using Gene Pix Pro 6.0 (Molecular Devices). Ratio based data normalization and selection of features were performed using Acuity 4.0 (Molecular Devices). All Spots with raw intensity less then 100U and less then twice the average background was ignored during normalization. Normalized data was filtered for the selection of features before further analysis. Only those spot were selected which contained a small percentage (<3) of saturated pixels, were not flagged bad or found absent (flags>0), had relatively uniform intensity and uniform background [Rgn R2 (635/532)>0.6] and were detectable above background (SNR>3). Analyzable spots in at least three of four biological replicates performed were retrieved for downstream analysis using Significant Analysis of Microarrays (SAM 3.0, Excel Add-In) [15], under the conditions of one class response and 100 permutations. Normalized log2 ratio (635/534) of four biological replicates with balanced dye-swaps was used for clustering. Microarrays were clustered using Acuity 4.0 (Molecular Devices). Gene symbols and FlyBase IDs against CG numbers were retrieved using FLIGHT (http://www.flight.licr.org/search/batch_genes.jsp) and GeneMerge (http://genemerge.bioteam.net/convertgenenames.html), respectively. GOTool Box [16] was used to retrieve overrepresented biological processes in upor downregulated genes, under the settings, ontology, biological process or molecular function; mode, all terms; reference, genome; evidence, all-all evidence; species, D. melanogaster; GO-stats; statistic test, hypergeometric (http://burgundy.cmmt.ubc.ca/GOToolBox/). Molecular functions overrepresented at hypergeometric distribution probability <0.05 after Bonferroni correction for multiple hypothesis testing were considered significant in the unbiased analysis. The protein interactome database BioGRID v 2.0 (http://www.thebiogrid.org/index.php) [17] was used in conjunction with the visualization software Osprey v. 1.2.0 (http://biodata.mshri.on.ca/osprey/servlet/Index) for examining overinteraction. The interactome database consisted of 7358 genes (vertices) and 24984 connections (edges). Overinteraction at hypergeometric distribution probability <0.05 after Bonferroni correction for multiple hypothesis testing was considered as significant. Enrichment for miRNA-target genes were analyzed using EMBL’s 2005 database (http://www.russell.embl-heidelberg.de/miRNAs/). As with protein interactome data, the software platform Osprey v. 1.2.0 was used for visualizing miRNA-target network. The miRNA-target database consisted of 4220 genes and 56 miRNAs. Target enrichment at hypergeometric distribution probability <0.05 after Bonferroni correction was considered as significant. For specifically examining enrichment of the GO biological process “axon guidance”, GOTool Box [16] was used as above except that Bonferroni correction was not applied. Human homologs of fly genes were identified using HomoloGene (http://www.ncbi.nlm.nih.gov/). Human homologs of axon guidance genes were depicted in the KEGG pathway for Homo sapiens (http://www.genome.jp/kegg/tool/color_pathway.html) because of unavailability of Drosophila pathway. Results and Discussion PTZ withdrawal induced behavioral alteration is ameliorated by AEDs Flies were treated with PTZ for 7 days, with AED for 3 days, and with drug free media for 4 days in chronological order, before climbing speed was measured. Flies treated with NF throughout served as controls. As reported earlier [11], PTZ withdrawal resulted in an enhanced speed (Figure 1). Notably, each of the five AEDs used ameliorated development of the PTZ withdrawal induced climbing abnormality (Figure 1). The climbing assay thus demonstrated that ETH, NaVP, LEV, GBP and VGB are all effective in post withdrawal model. This was in contrast to previously described chronic PTZ model in which only NaVP and LEV, not ETH, GBP and VGB, were found to be effective in ameliorating reduction in climbing speed caused by the convulsant drug [11]. In rodent epilepsy model, PTZ kindling is induced by repeated injection of the GABA antagonist over several weeks, resulting in partially and fully kindled animals; once kindled, the state of behavioral hyperexcitability persists for up to several weeks after discontinuation of PTZ treatment [18, 19]. The pharmacology of kindled and postkindled rodent epilepsy models is known to differ [20]. It was thus notable that the fly model showed different behavioral pharmacology between chronic PTZ treatment and post-PTZ withdrawal regimes. PTZ withdrawal induced transcriptomic alteration We generated expression profiles of fly heads at three time points – 1, 3 and 7 day after PTZ withdrawal, i.e., on 8, 10 and 14 day from the start of PTZ treatment. These time-points represented beginning, latent, and behaviorally expressive phase of post-PTZ withdrawal, in that order. The full microarray data set has been deposited in the Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/) under accession series GSE7156. Clustering of microarrays showed the beginning phase behaving as an outgroup (Figure 2). A dynamic change in the transcriptome was thus found to underlie the fly model. We next identified differentially expressed genes. Previously, control microarray experiment comparing NF versus NF treated flies using the method used here was reported to detect no differentially expressed gene below 96% False Discovery Rate (FDR) [11]. In the present microarrays, we used a cut-off of 15% FDR (for gene lists, see additional file 1). In a time-series, differentially expressed genes in adjacent timepoints would be expected to match significantly in a direction-specific manner. We examined this using hypergeometric distribution probability of overlap, assuming population sizes of 10000, approximately the number of unique genes in the arrays. Whereas 8 and 10 day genes did not match significantly, a direction specific match was indeed observed between 10 and 14 day. The upand down-regulated genes in these later time-points showed significant overlap, with p values of 1.72E-09 and 3.23E11, respectively. Further analysis showed enrichment of several GO molecular functions in differentially expressed genes at these two time-points (Table 1). Overall, post-PTZ withdrawal regime was characterized by upregulation of genes overrepresenting transcription regulatory (Figure 3) and GTPase regulator activities, and downregulation of genes enriched in protein binding, oxidoreductase and electron carrier activities. It was interesting to note here the already existing gene expression evidence for upregulated Gprotein signaling and transcriptional regulation activities during post-epileptogenesis in rodent models [21-24]. Amelioration of PTZ withdrawal induced transcriptomic alteration by AEDs Since AEDs ameliorated development of locomotor abnormality post-PTZ withdrawal (Figure 1), it was of interest to examine if AEDs act through G-protein signaling and transcriptional regulation in our fly model. For this, we generated microarray expression profiles of flies treated with AEDs post-PTZ withdrawal. Flies were treated with PTZ for seven days, following which they were first treated for 3 days with ETH, GBP, VGB, NaVP or LEV and then with drug free media for next four days. Microarray expression profiles of heads were generated from flies collected on 10 and 14 day from the beginning of PTZ treatment. Clustering of expression profiles was found to be consistent with the time-points (Figure 4). This suggested similarity in AEDs’ action. The full microarray data set has been deposited in the Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/) under accession series GSE10984 (ETH), GSE10985 (GBP), GSE10986 (NaVP), GSE10987 (VGB) and GSE10988 (LEV). Differentially expressed genes in AED expression profiles were retrieved at 15% FDR cut-off (for gene lists, see additional file 2). Genes were predominantly downregulated by AEDs. Further analysis showed enrichment of various GO molecular functions in the differentially expressed genes (Table 2). Whereas these functions included oxidoreductase and electron carrier activities found enriched in AED-untreated group (Table 1), it was striking that unlike the latter none of the AED-treated group showed overrepresentation of transcription regulatory and GTPase regulator activities. Altogether, our transcriptomic analysis suggested that post-PTZ withdrawal induced alterations in transcription regulatory and GTPase regulator activities are ameliorated by AEDs. Interestingly, microarray expression profiling of rat brain has previously shown that NaVP regulates transcription regulatory and G-protein signaling [25]. Cumulatively, our results suggested that post-PTZ withdrawal in Drosophila models post-epileptogenesis. Protein interactome and miRNA-target based analyses of post-PTZ withdrawal genes We next used the available protein interactome map of Drosophila [17] to examine if it can provide directionality to gene expression changes (additional file 1) post-PTZ withdrawal. For this, we first analyze differentially expressed 10 day genes under the GO categories transcription regulator activity and GTPase regulator activity (Table 1) for overinteraction. Out of total 78 genes under these categories, 31 showed one or more interacting partners within the entire set of differentially expressed post-PTZ withdrawal time-series genes. Considering the total number of genes in the database as 7358, two of the 31 genes tested, the GTPase activating protein RanGap (CG9999) and the transcription factor ovo (CG6824), showed significant overinteraction. Interacting with each other, these genes showed three and four partners, respectively. Considering six and 11 total interacting partners in the database, these interactions were significant (hypergeometric distribution p=0.046 and 0.027, Bonferroni adjusted, in that order). Considering this finding as important, we segregated all the seven unique interacting partners in the RanGap and ovo sub-network on post-PTZ withdrawal time scale. These seven genes spread to all the three time-points (Figure 5). They are highly conserved genes. For example, all except ovo has a human homolog KRAS (Ras85D), RANGAP1 (RanGap), DHX9 (mle), TCF7L2 (pan), PELI2 (Pli), and LMO4 (CG5708). The beginning time-point mammalian homolog KRAS maps to various pathways implicated in epileptogenesis axon guidance, long-term depression, and MAPK signaling [26-28]. Although the rest of the five homologs did not map to these pathways, it is notable that the last time-point homolog LMO4 has recently been found to interact with a small GTPase activating repulsive guidance molecule [29]. Also, downregulation of LMO4 has been shown to neutralize repulsive activity of the molecule [29]. To test the hypothesis that axon guidance abnormality underlies our post-PTZ withdrawal model, we examined GO biological process enrichment in differentially expressed genes at all three timepoints combined together (additional file 1). Remarkably, axon guidance was one among several other processes such as neurogenesis, axonogenesis etc., that showed enrichment (for process list, see additional file 3). Both upand down-regulation of axon guidance genes was observed (Table 3). Notably, many of the processes enriched in our model including axon guidance are known to be associated with rodent models of kindling epilepsy and status epilepticus [30]. In the axon guidance pathway, membrane receptors predominantly represented differentially expressed genes (Figure 6). Notably, axon guidance receptors have been implicated in rodent epilepsy models [31, 32]. Overall, protein-interaction analysis provided evidence for an altered axon guidance pathway in the post-PTZ withdrawal model. Most of the axon guidance genes belonged to 10 day time-point (Table 4). All of them were upregulated. To further test the axon guidance hypothesis, we used available miRNA-target database and examined if 10 day upregulated axon guidance genes are enriched for miRNA-targets. A total of 38 miRNA showed one or more targets in the gene set. Significant enrichment of targets was observed for mir-34 (Figure 7). Of the total 14 10 day genes, five were the targets of mir-34. Considering 178 total mir-34 targets in the database containing 4220 genes, this enrichment was significant (hypergeometric distribution p=0.006, Bonferroni adjusted). The miRNA mir-34 has previously been shown to be regulated by the steroid hormone ecdysone and the terpenoid juvenile in Drosophila [33]. Hormone-dependent rewiring of neurons is a conserved feature of brain development and plasticity [34]. Hormones are known to regulate adult behavior and associated neuronal plasticity in insects [34]. In human epilepsy also, steroid hormones are considered as pharmacoactive compounds influencing seizure threshold [35]. Given this, it was interesting to find mir-34 target enrichment in our model. Cumulatively, miRNA-target analysis supported a role of neuronal remodeling in post-PTZ withdrawal model.